Story pitch:

I’m analyzing the percentage change in childhood vaccination rates in Michigan by county across the past decade, with a focus on measles as Michigan continues to see more cases of the disease.

<script src="https://public.flourish.studio/resources/embed.js"></script>
<noscript><img src="https://public.flourish.studio/visualisation/23592389/thumbnail" width="100%" alt="map visualization" /></noscript>
<script src="https://public.flourish.studio/resources/embed.js"></script>
<noscript><img src="https://public.flourish.studio/visualisation/23636168/thumbnail" width="100%" alt="map visualization" /></noscript>

After pulling state-level data for the measles, mumps and rubella vaccine directly from the Michigan Department of Health and Human Services website, I found that the infant vaccination rate in Michigan dropped slightly in the past five years by about six percent, while the teenage MMR rate has stayed relatively steady in the past year, which was all the data that was available on the site. Data was not available prior to 2020 for individual antigens, such as the MMR vaccine.

Vaccine waivers

Based on my research, there are two main thresholds for children to receive their recommended vaccinations: Kindergarten and seventh grade. Parents and families are able to submit vaccine waivers due to religion, medical reasons and philosophy/personal reasons. After looking at 2023 data on vaccine waivers, I found that Houghton County had high rates of waivers at both grade thresholds. This is likely due to the small number of students attending school in both of these grades. I’d like to perform further analysis on the types of vaccine waivers and how that differs by county, but was not able to get to this in time.

Median income and vaccine rates

I visualized income and childhood and adolescent vaccine rates, but based on the correlation coefficient analysis, it seems like these two do not have significant correlation.

I performed a correlation coefficient analysis on both the percent change in election results (percentage point difference between parties) and the percentage change in the vaccination rates for each of the two series for children and teenagers, respectively, and found that there was little correlation. I also performed the same analysis on aggregated vaccine series data and election results for each election year and found that there was also no significant correlation here. Because of this, I’m moving away from the angle of looking at election results to see if counties that vote more Republican line up with a drop in vaccination rates across counties. Instead, I’m shifting my angle to just look at overall changes in vaccination based on the vaccination report card data, to determine the counties that have dropped the most in their recommended vaccination coverage for children and teenagers.

I initially tried to aggregate the four vaccine series, but I do not think that is the best method to go about this analysis, as they cover different age ranges. Instead, I’m aggregating the vaccine rates based on the age range that is recommended to get them — two of the series are specifically for children aged 19 months to 35 months, and the two others are recommended for adolescents aged 13 to 17.

When analyzing each of the four childhood series individually, I found that two counties in Michigan — Oscoda and Keweenaw — have seen the highest drops in immunization rates for three of the four major vaccine series, and consistently fell in the bottom three for lowest MMR vaccination rates in the state. These two counties also have some of the lowest populations in the state, which could potentially indicate that more rural areas are seeing larger drops in their vaccination rates for children.

I originally joined the vaccination report card data to Census ACS data from 2023 on population for each county, and am working on joining election data to this to see if there’s a pattern between vaccination rates and how people voted in the 2024 and 2016 elections. I’d like to see the proportion of votes allocated toward the Republican Party and the Democratic Party for each of the counties in Michigan, and I want to perform a more robust analysis for these counties beyond what I found of the two counties that immediately stand out. I also plan on looking at the vaccination rates for the MMR vaccine alone. I’ve found the csv for this from the Michigan government and will load it in as I continue work on this story. The measles vaccine, which is included in each of the four vaccine series I’m analyzing, is primarily administered to young children and adolescents. There have been 9 reported measles cases this year in Michigan so far, in Ingham, Kent, Macomb, Marquette, Montcalm and Oakland Counties, making this the worst year in decades for measles in Michigan.

Load software libraries 

#install.packages("tidyverse")
#install.packages("janitor")
#install.packages("readxl")
#install.packages("tidycensus")
#install.packages("DT")
library(tidyverse)
library(janitor)
library(readxl)
library(tidycensus)
library(DT)
library(sf)
library(gridExtra)
library(zoo)
library(ellmer)
library(showtext)

Load data + initial cleaning

# The report card data provides estimates for immunization coverage at the county level for vaccines recommended for children (19-35 months) and adolescent (13-17-year-olds).
# Load vaccine report card data.
report_card_data <- read.csv("../data/Report_card_All_data_20250313.csv") |> 
    clean_names()

# Load geoJSON for mapping.
mi_counties <- read_sf("../data/County.geojson") %>% clean_names()


# How many counties?
report_card_data$cntyname %>% 
  unique()
##  [1] "Alcona"         "Alger"          "Allegan"        "Alpena"        
##  [5] "Antrim"         "Arenac"         "Baraga"         "Barry"         
##  [9] "Bay"            "Benzie"         "Berrien"        "Branch"        
## [13] "Calhoun"        "Cass"           "Charlevoix"     "Cheboygan"     
## [17] "Chippewa"       "Clare"          "Clinton"        "Crawford"      
## [21] "Delta"          "Detroit"        "Dickinson"      "Eaton"         
## [25] "Emmet"          "Genesee"        "Gladwin"        "Gogebic"       
## [29] "Grand Traverse" "Gratiot"        "Hillsdale"      "Houghton"      
## [33] "Huron"          "Ingham"         "Ionia"          "Iosco"         
## [37] "Iron"           "Isabella"       "Jackson"        "Kalamazoo"     
## [41] "Kalkaska"       "Kent"           "Keweenaw"       "Lake"          
## [45] "Lapeer"         "Leelanau"       "Lenawee"        "Livingston"    
## [49] "Luce"           "Mackinac"       "Macomb"         "Manistee"      
## [53] "Marquette"      "Mason"          "Mecosta"        "Menominee"     
## [57] "Michigan"       "Midland"        "Missaukee"      "Monroe"        
## [61] "Montcalm"       "Montmorency"    "Muskegon"       "Newaygo"       
## [65] "Oakland"        "Oceana"         "Ogemaw"         "Ontonagon"     
## [69] "Osceola"        "Oscoda"         "Otsego"         "Ottawa"        
## [73] "Presque Isle"   "Roscommon"      "Saginaw"        "Sanilac"       
## [77] "Schoolcraft"    "Shiawassee"     "St. Clair"      "St. Joseph"    
## [81] "Tuscola"        "Van Buren"      "Washtenaw"      "Wayne"         
## [85] "Wexford"
# Michigan has 83 counties. There are 85 unique county values, because this data set includes Detroit City and the State of Michigan's vaccination report cards.
# I'm creating a column based on an index so I can plot this and analyze by time. 
report_card_data <- report_card_data %>% 
  mutate(
     quarter_num = as.numeric(sub("Q", "", quarter)),
     quarter_index = year + (quarter_num - 1) / 4,
    quarter_date = make_date(year, (quarter_num - 1) * 3 + 1, 1)
  )

Immunization waivers by school: kindergarten

school_vaccine_k <- read_excel("../data/kindergarten-immunization-school.xlsx", skip = 7) %>% clean_names()
## New names:
## • `n` -> `n...10`
## • `%` -> `%...11`
## • `n` -> `n...12`
## • `%` -> `%...13`
## • `n` -> `n...14`
## • `%` -> `%...15`
## • `n` -> `n...16`
## • `%` -> `%...17`
school_vaccine_k %>% 
  filter(type == "Pub") %>% 
  group_by(county) %>% 
  summarise(weighted_pct = sum(percent_comp * n) / sum(n)) %>% 
  arrange(weighted_pct)
## # A tibble: 82 × 2
##    county      weighted_pct
##    <chr>              <dbl>
##  1 Antrim              81.7
##  2 Oscoda              83.6
##  3 Manistee            83.8
##  4 Lapeer              84.4
##  5 Houghton            84.5
##  6 Wexford             85.6
##  7 Emmet               85.8
##  8 Detroit             85.8
##  9 Montmorency         86.0
## 10 Newaygo             86.0
## # ℹ 72 more rows
# percent of total vaccine waivers by county
school_vaccine_k %>% 
  group_by(county) %>% 
  summarize(total_waiver_pct = sum(percent_11 * n_10)/ sum(n_10)) %>% 
  arrange(desc(total_waiver_pct))
## # A tibble: 82 × 2
##    county      total_waiver_pct
##    <chr>                  <dbl>
##  1 Benzie                  44.4
##  2 Ogemaw                  27.8
##  3 Ontonagon               25  
##  4 Houghton                23.6
##  5 Alpena                  21.7
##  6 St. Joseph              21.2
##  7 Montmorency             21.1
##  8 Antrim                  19.6
##  9 Lapeer                  19.4
## 10 Kalamazoo               18.4
## # ℹ 72 more rows
waiver_rate_k <- school_vaccine_k %>%
  filter(type == "Pub") %>% 
  group_by(county) %>%
  summarize(
    total_waivers = sum(n_10, na.rm = TRUE),
    total_students = sum(n, na.rm = TRUE),
    waiver_rate = round(100 * total_waivers / total_students, 1)
  ) %>%
  arrange(desc(waiver_rate))

write_csv(waiver_rate_k, "../data/waiver-rate_k.csv")

datatable(waiver_rate_k, options = list(pageLength = 10))
#install.packages("showtext")
font_add_google("Roboto Condensed", "roboto")
showtext_auto()


mi_counties %>% 
  inner_join(waiver_rate_k, by = c("name" = "county")) %>% 
  ggplot() +
  labs(
    title = "Houghton, Antrim and Emmet Counties had the highest rates of
kindergarten immunization waivers in Michigan, each over 13%",
  subtitle = "None of Luce County's 50 kindergarten-aged students filed vaccine waivers."
  )+
  geom_sf(aes(fill = waiver_rate)) +
  scale_fill_gradient(
    low = "#efefd0",
    high = "#ff6b35"
  ) +
  theme(
    text = element_text(family = "roboto"),
    plot.title = element_text(size = 14, face = "bold")
  )

Immunization waivers by school: seventh grade

school_vaccine_seventh <- read_excel("../data/seventh-grade-waivers.xlsx", skip = 7) %>% clean_names()
## New names:
## • `n` -> `n...10`
## • `%` -> `%...11`
## • `n` -> `n...12`
## • `%` -> `%...13`
## • `n` -> `n...14`
## • `%` -> `%...15`
## • `n` -> `n...16`
## • `%` -> `%...17`
school_vaccine_seventh %>% 
  filter(type == "Pub") %>% 
  group_by(county) %>% 
  summarise(weighted_pct = sum(percent_comp * n) / sum(n)) %>% 
  arrange(weighted_pct)
## # A tibble: 84 × 2
##    county                weighted_pct
##    <chr>                        <dbl>
##  1 No County Affiliation         66.7
##  2 Oscoda                        78.5
##  3 Houghton                      80.3
##  4 Alcona                        81.1
##  5 Leelanau                      82.2
##  6 Montmorency                   83.7
##  7 Manistee                      83.8
##  8 Antrim                        84.3
##  9 Arenac                        84.4
## 10 Lapeer                        84.5
## # ℹ 74 more rows
# percent of total vaccine waivers by county
school_vaccine_seventh %>% 
  group_by(county) %>% 
  summarize(total_waiver_pct = sum(percent_11 * n_10)/ sum(n_10)) %>% 
  arrange(desc(total_waiver_pct))
## # A tibble: 84 × 2
##    county                total_waiver_pct
##    <chr>                            <dbl>
##  1 No County Affiliation             33.3
##  2 Leelanau                          27.4
##  3 Houghton                          26.2
##  4 Oscoda                            24.9
##  5 Montcalm                          21.4
##  6 Alpena                            20.6
##  7 Clare                             17.6
##  8 Eaton                             16.4
##  9 Montmorency                       16.0
## 10 St. Joseph                        15.9
## # ℹ 74 more rows
waiver_rate_7th <- school_vaccine_seventh %>%
  filter(type == "Pub") %>% 
  group_by(county) %>%
  summarize(
    total_waivers = sum(n_10, na.rm = TRUE),
    total_students = sum(n, na.rm = TRUE),
    waiver_rate = round(100 * total_waivers / total_students, 1)
  ) %>%
  arrange(desc(waiver_rate))

write_csv(waiver_rate_7th, "../data/waiver_rate_7th.csv")

datatable(waiver_rate_7th, options = list(pageLength = 10))
#install.packages("showtext")
font_add_google("Roboto Condensed", "roboto")
showtext_auto()


mi_counties %>% 
  inner_join(waiver_rate_7th, by = c("name" = "county")) %>% 
  ggplot() +
  labs(
    title = "Oscoda, Houghton, Montmorency and Alger Counties had 
the highest rates of seventh grade immunization waivers 
in Michigan",
  subtitle = "About three percent of seventh grade students in Detroit City filed vaccination 
waivers in 2023."
  )+
  geom_sf(aes(fill = waiver_rate)) +
  scale_fill_gradient(
    low = "#efefd0",
    high = "#004e89"
  ) +
  theme(
    text = element_text(family = "roboto"),
    plot.title = element_text(size = 14, face = "bold")
  )

Percent change analysis for four childhood vaccination series over the past decade

### PERCENT CHANGE IN VACCINATION COVERAGE

# For this analysis, I'm focusing on the vaccine series 4313314, 43133142, 132321, 1323213, and the MMR vaccine. These are the series Michigan recommends for children and adolescents — the first two are for babies from 19-35 months, and the latter two are for teens.

4313314 series:

  • 4 or more doses of DTaP/DTP/DT
  • 3 or more doses of Polio
  • 1 or more dose of MMR
  • 3 or more doses of Hib
  • 3 or more doses of HepB
  • 1 or more dose of Varicella
  • 4 or more doses of PCV.
# Let's look at the 4313314 series: "4 or more doses of DTaP/DTP/DT, 3 or more doses of Polio, 1 or more dose of MMR, 3 or more doses of Hib, 3 ore more doses of HepB, 1 or more dose of Varicella, 4 or more doses of PCV." 
# Coverage for the 4313314 vaccine decreased for all but three counties in Michigan. Oscoda and Keweenaw Counties nearly halved their vaccine coverage for this series. 

series_1 <- report_card_data %>% 
  group_by(cntyname) %>% 
  filter(quarter_index %in% c(2014.75, 2024.75)) %>% 
  select(cntyname, quarter_index, covg4313314) %>% 
  pivot_wider(
    names_from = quarter_index,
    values_from = covg4313314,
    names_prefix = "Q_"
  ) %>% 
  mutate(pct_change_2014_2024 = 100 * (Q_2024.75 - Q_2014.75) / Q_2014.75) %>% 
  arrange(pct_change_2014_2024)

datatable(series_1, options = list(pageLength = 10))

43133142 series:

  • 4 or more doses of DTaP/DTP/DT
  • 3 or more doses of Polio
  • 1 or more dose of MMR
  • 3 or more doses of Hib
  • 3 or more doses of Hep B
  • 1 or more dose of Varicella
  • 4 or more doses of PCV
  • 2 or more doses of HepA
## Again, Oscoda and Keweenaw are in the top two counties that have nearly halved their vaccine coverage for this series.
series_2 <- report_card_data %>% 
  group_by(cntyname) %>% 
  filter(quarter_index %in% c(2014.75, 2024.75)) %>% 
  select(cntyname, quarter_index, covg43133142) %>% 
  pivot_wider(
    names_from = quarter_index,
    values_from = covg43133142,
    names_prefix = "Q_"
  ) %>% 
  mutate(pct_change_2014_2024 = 100 * (Q_2024.75 - Q_2014.75) / Q_2014.75) %>% 
  arrange(pct_change_2014_2024)

datatable(series_2, options = list(pageLength = 10))

132321 series:

  • 1 or more doses of Tdap
  • 3 or more doses of Polio
  • 2 or more doses of MMR
  • 3 or more doses of HepB
  • 2 or more doses of Varicella
  • 1 or more dose of MenACWY
# This series did not have as high percent change decreases, but Keewenaw is in the top five here too.
series_3 <- report_card_data %>% 
  group_by(cntyname) %>% 
  filter(quarter_index %in% c(2014.75, 2024.75)) %>% 
  select(cntyname, quarter_index, covg132321) %>% 
  pivot_wider(
    names_from = quarter_index,
    values_from = covg132321,
    names_prefix = "Q_"
  ) %>% 
  mutate(pct_change_2014_2024 = 100 * (Q_2024.75 - Q_2014.75) / Q_2014.75) %>% 
  arrange(pct_change_2014_2024)

datatable(series_3, options = list(pageLength = 10))

1323213 series:

  • 1 or more doses of Tdap
  • 3 or more doses of Polio
  • 2 or more doses of MMR
  • 3 or more doses of HepB
  • 2 or more doses of Varicella
  • 1 or more dose of MenACWY
  • HPV complete (with 2 or 3 doses (Males & Females))
# Nearly all counties have increased their coverage of this series. I'm wondering the reason behind this, and if it is significant. Could this change be because this is a new vaccine series that only started being implemented recently? Would love to talk to a Michigan epidemiologist who works for the state DHHS or some other expert to understand this. 
series_4 <- report_card_data %>% 
  group_by(cntyname) %>% 
  filter(quarter_index %in% c(2014.75, 2024.75)) %>% 
  select(cntyname, quarter_index, covg1323213) %>% 
  pivot_wider(
    names_from = quarter_index,
    values_from = covg1323213,
    names_prefix = "Q_"
  ) %>% 
  mutate(pct_change_2014_2024 = 100 * (Q_2024.75 - Q_2014.75) / Q_2014.75) %>% 
  arrange(pct_change_2014_2024)

datatable(series_4, options = list(pageLength = 10))

Analysis by age of series

The 431XXXX series are administered to children 19-35 months, while the 132XXXX series are for adolescents aged 13-17 years. Let’s do an unweighted mean of vaccination rates by the age range of the vaccine series.

# Aggregating vaccine rates for just the infant series.
baby_vax <- report_card_data %>% 
  mutate(vax_agg = ((covg4313314 + covg43133142)/2)) %>% 
  select(cntyname, year, vax_agg) %>% 
  group_by(cntyname, year) %>% 
  summarize(mean_vax_rate = mean(vax_agg, na.rm = TRUE)) 

# Infant vaccine series rate over past decade.
baby_vax_change <- baby_vax %>% 
  pivot_wider(
    names_from = year,
    values_from = mean_vax_rate,
    names_prefix = "Y_"
  ) %>% 
  mutate(
    pct_change_14_24 = ((Y_2024 - Y_2014)/Y_2014) * 100
  ) %>% 
  select(cntyname, pct_change_14_24) %>% 
  arrange(pct_change_14_24)


write_csv(baby_vax_change, "../data/infant-vax-pct-change.csv")

# Mapping percentage change of vaccination rates for infants.
# mi_counties %>% 
#   inner_join(baby_vax_change, by = c("name" = "cntyname")) %>% 
#   ggplot() +
#   geom_sf(aes(fill = pct_change_14_24)) +
#   scale_fill_gradient(
#     high = "#FFFFC5",
#     low = "#e51d37") +
#   labs(title = "Rates of recommended vaccines for children aged 19 to 35 months
# decreased by 46% in Keweenaw County during the past decade,
# the highest in Michigan") +
#   theme_void()

# Aggregating vaccine rates for just the adolescent series.
teen_vax <- report_card_data %>% 
  mutate(vax_agg = ((covg132321 + covg1323213)/2)) %>% 
  select(cntyname, year, vax_agg) %>% 
  group_by(cntyname, year) %>% 
  summarize(mean_vax_rate = mean(vax_agg, na.rm = TRUE)) 

# Adolescent vaccine rate over past decade.
# The vaccine rate for the teen series has increased in nearly every county in the past ten years, primarily because of the 1323213 series, which saw an increase in nearly every county as well. I wonder if this could be because it's a newer series that only recently began to be implemented, but I need to talk to an expert about this.
teen_vax_change <- teen_vax %>% 
  pivot_wider(
    names_from = year,
    values_from = mean_vax_rate,
    names_prefix = "Y_"
  ) %>% 
  mutate(
    pct_change_14_24 = ((Y_2024 - Y_2014)/Y_2014) * 100
  ) %>% 
  select(cntyname, pct_change_14_24) %>% 
  arrange(pct_change_14_24)

write_csv(teen_vax_change, "../data/teen-vax-pct-change.csv")

# Mapping percentage change of teenage vaccination series from 2014 to 2024.
# mi_counties %>% 
#   inner_join(teen_vax_change, by = c("name" = "cntyname")) %>% 
#   ggplot() +
#   geom_sf(aes(fill = pct_change_14_24)) +
#   theme_void()

Analyze percentage of waiver types: philosophy (personal reasons), medical, religion.

# get info from waiver data table

table_image <- content_image_file("../data/table-ss.png")

chat <- chat_groq(
  model = "meta-llama/llama-4-scout-17b-16e-instruct"
)

# waiver_data <- chat$chat("turn the data from this image into a structured csv format", table_image)

MMR vaccine data: 19-35 months.

The MMR vaccine’s coverage rate in young children has decreased by 6.8% since January 2020. 
infant_vax_2020_25 <- read_excel("../data/mcir-vax-data/Profile-by-County-19-35-Months-2020-Present.xlsx", sheet = 2) %>% clean_names

# Filtering the statewide MMR vaccine rates for Michigan.
infant_mmr_2020_25 <- infant_vax_2020_25 %>% 
  filter(antigen == "1+ MMR") %>% 
  pivot_longer(
    cols = jan2020:apr2025,
    names_to = "month",
    values_to = "vax_rate"
  ) %>% 
  mutate(month_year = as.yearmon(month, format = "%b%Y")) %>% 
  select(-month) 

# Plotting change over time.
infant_mmr_2020_25 %>% 
  ggplot() +
  geom_line(aes(x=month_year, y = vax_rate)) +
  expand_limits(y=0)

# Percent change from January 2020 to April 2025.
100 * ((79.5-85.3)/85.3)
## [1] -6.799531

MMR vaccine data: 13-17 years. (Data is only available for May 2024 to April 2025)

MMR vaccine rates for teenagers has stayed relatively the same in the past year, increasing by 0.23% in that time. 

# Filtering just statewide MMR vaccine data for adolescents.
teen_mmr_2020_25 <- read_excel("../data/mcir-vax-data/Profile-by-County-13-17-Years-2020-Present.xlsx", sheet = 2) %>% clean_names %>% 
  filter(antigen == "2+ MMR") %>% 
  pivot_longer(
    cols = may2024:apr2025,
    names_to = "month",
    values_to = "vax_rate"
  ) %>% 
  mutate(month_year = as.yearmon(month, format = "%b%Y")) %>% 
  select(-month) 

100 * ((86.5-86.3)/86.3)
## [1] 0.2317497

Median income and vaccine rates: childhood series.

mi_income <- get_acs(geography = "county", 
              variables = c(medincome = "B19013_001"), 
              state = "MI", 
              year = 2023) %>% clean_names() %>% arrange(estimate)
## Getting data from the 2019-2023 5-year ACS
mi_income$name <- sub(" County.*", "", mi_income$name)

baby_vax_and_income <- baby_vax_change %>% 
  inner_join(mi_income, by=c("cntyname" = "name")) %>% 
  select(cntyname, pct_change_14_24, estimate)

write.csv(baby_vax_and_income, "../data/infant-vax-and-income.csv")

cor(baby_vax_and_income$pct_change_14_24, baby_vax_and_income$estimate)
## [1] 0.2309638
# No clear correlation between childhood vaccination rate and median income.

Median income and vaccine rates: adolescent series.

teen_vax_and_income <- teen_vax_change %>% 
  inner_join(mi_income, by=c("cntyname" = "name")) %>% 
  select(cntyname, pct_change_14_24, estimate) %>% 
  arrange(estimate)

write.csv(teen_vax_and_income, "../data/teen-vax-and-income.csv")

cor(teen_vax_and_income$pct_change_14_24, teen_vax_and_income$estimate)
## [1] 0.2911488
# No clear correlation between teen vax rate and median income.

Population analysis

# Let's look at populations for these counties.
# Load 2023 ACS population data from tidycensus.
vars <- load_variables(2023, "acs1") 

# 2023 Census population by county.
pop_by_county <- get_acs(
  geography = "county",
  variables = "B01001A_001",
  state = "MI",
) %>% clean_names()

# Cleaning census data set's names to get rid of extra text like "County."
pop_by_county$name <- sub(" County.*", "", pop_by_county$name)

report_card_with_pop <- report_card_data %>% 
  inner_join(pop_by_county, join_by(cntyname == name)) %>% 
  select(
    cntyname, estimate, quarter_index, covg4313314, covg43133142, covg132321, covg1323213 
  )

datatable(pop_by_county, options = list(pageLength = 10))
# Keewenaw is the county with the smallest population, with 2,021 people. 
report_card_with_pop %>%
  filter(quarter_index == "2023.75") %>% 
  group_by(cntyname) %>% 
  arrange(estimate)
## # A tibble: 83 × 7
## # Groups:   cntyname [83]
##    cntyname    estimate quarter_index covg4313314 covg43133142 covg132321
##    <chr>          <dbl>         <dbl>       <dbl>        <dbl>      <dbl>
##  1 Keweenaw        2021         2024.        43.8         37.5       67.7
##  2 Luce            4909         2024.        72.1         50         76.5
##  3 Ontonagon       5487         2024.        68.9         66.7       76.4
##  4 Baraga          5991         2024.        63.6         54.5       77.6
##  5 Schoolcraft     6820         2024.        76.8         65.3       78.8
##  6 Alger           7234         2024.        67           59.1       71.9
##  7 Mackinac        7827         2024.        61           44.8       73.4
##  8 Oscoda          7856         2024.        43           29.9       59.5
##  9 Montmorency     8802         2024.        64.3         46         76.8
## 10 Alcona          9727         2024.        69.6         54.9       73.5
## # ℹ 73 more rows
## # ℹ 1 more variable: covg1323213 <dbl>
# report_card_data %>% 
#   filter(cntyname == "Keweenaw") %>% 
#   ggplot(aes(x=quarter_index)) +
#   geom_line(aes(y = covg4313314, color = "covg4313314")) +
#   geom_line(aes(y = covg43133142, color = "covg43133142")) +
#   geom_line(aes(y = covg132321, color = "covg132321")) +
#   geom_line(aes(y = covg1323213, color = "covg1323213")) +
#   labs(title = "Coverage Over Time in Keweenaw County",
#        y = "Coverage Value",
#        x = "Quarter",
#        color = "Coverage Type") +
#   theme_minimal()

Aggregate election analysis

election_16 <- read.csv("../data/2016_US_County_Level_Presidential_Results.csv") %>%  clean_names()
election_20 <- read.csv("../data/2020_US_County_Level_Presidential_Results.csv") %>% clean_names()
election_24 <- read.csv("../data/2024_US_County_Level_Presidential_Results.csv") %>% clean_names()

# Filtering for just Michigan counties
mi_election_16 <- election_16 %>% 
  filter(state_abbr == "MI")
mi_election_20 <- election_20 %>% 
  filter(state_name == "Michigan")
mi_election_24 <- election_24 %>% 
  filter(state_name == "Michigan")

# cleaning the names so they don't have "County" appended, since we're joining based on name
mi_election_16$county_name <- sub(" County.*", "", mi_election_16$county_name)
mi_election_20$county_name <- sub(" County.*", "", mi_election_20$county_name)
mi_election_24$county_name <- sub(" County.*", "", mi_election_24$county_name)

# Let's look at the average vaccine coverage by county in just election years for INFANT vaccines.
 vaccine_election_2016 <- baby_vax %>%
   filter(year == 2016) %>%
   group_by(cntyname) %>%
   pivot_wider(
     names_from = year,
     values_from = mean_vax_rate,
     names_prefix = "Y_"
   ) %>%
   inner_join(mi_election_16, by = c("cntyname" = "county_name")) %>%
   select(cntyname, Y_2016, total_votes, per_dem, per_gop) %>%
   mutate(
     per_point_diff = (per_gop - per_dem)
   )

vaccine_election_2020 <- baby_vax %>%
  filter(year == 2020) %>%
  group_by(cntyname) %>%
  pivot_wider(
    names_from = year,
    values_from = mean_vax_rate,
    names_prefix = "Y_"
  ) %>%
   inner_join(mi_election_20, by = c("cntyname" = "county_name")) %>%
   select(cntyname, Y_2020, total_votes, per_dem, per_gop, per_point_diff)

 vaccine_election_2024 <- baby_vax %>%
   filter(year == 2024) %>%
   group_by(cntyname) %>%
   pivot_wider(
     names_from = year,
     values_from = mean_vax_rate,
     names_prefix = "Y_"
   ) %>%
   inner_join(mi_election_24, by = c("cntyname" = "county_name")) %>%
   select(cntyname, Y_2024, total_votes, per_dem, per_gop, per_point_diff)

Mapping election and vaccine results for 2024

# Plotting elections + vaccine results using a geoJSON of Michigan counties
# install.packages("gridExtra")

# joining geoJSON to aggregate childhood vaccination data 
vaccines_24 <- mi_counties %>% 
  inner_join(vaccine_election_2024, by = c("name" = "cntyname")) %>% 
  ggplot() +
  geom_sf(aes(fill = Y_2024)) +
  theme_void()

mi_election_map_24 <- mi_counties %>% 
  inner_join(vaccine_election_2024, by = c("name" = "cntyname")) %>% 
  ggplot() +
  geom_sf(aes(fill = per_point_diff)) +
  scale_fill_gradient(
    low = "#0041C2",
    high = "#e51d37"
  ) +
  theme_void()

grid.arrange(vaccines_24, mi_election_map_24, ncol = 2)

Election and vaccine rates in 2020

vaccines_20 <- mi_counties %>% 
  inner_join(vaccine_election_2020, by = c("name" = "cntyname")) %>% 
  ggplot() +
  geom_sf(aes(fill = Y_2020)) +
  theme_void()

mi_election_map_20 <- mi_counties %>% 
  inner_join(vaccine_election_2020, by = c("name" = "cntyname")) %>% 
  ggplot() +
  geom_sf(aes(fill = per_point_diff)) +
  scale_fill_gradient(
    low = "#0041C2",
    high = "#e51d37"
  ) +
  theme_void()

grid.arrange(vaccines_20, mi_election_map_20, ncol = 2)

Election and vaccine rates in 2016

vaccines_16 <- mi_counties %>% 
  inner_join(vaccine_election_2016, by = c("name" = "cntyname")) %>% 
  ggplot() +
  geom_sf(aes(fill = Y_2016)) +
  theme_void()

mi_election_map_16 <- mi_counties %>% 
  inner_join(vaccine_election_2016, by = c("name" = "cntyname")) %>% 
  ggplot() +
  geom_sf(aes(fill = per_point_diff)) +
  scale_fill_gradient(
    low = "#0041C2",
    high = "#e51d37"
  ) +
  theme_void()

grid.arrange(vaccines_16, mi_election_map_16, ncol = 2)

Correlation coefficients

# Correlation coefficients between the vaccine rates and the percentage point difference in that year's election. (GOP% - DEM%)
cor(vaccine_election_2020$Y_2020, vaccine_election_2020$per_point_diff)
## [1] -0.2696495
cor(vaccine_election_2024$Y_2024, vaccine_election_2024$per_point_diff)
## [1] -0.3843817
cor(vaccine_election_2016$Y_2016, vaccine_election_2016$per_point_diff)
## [1] -0.2238821
# No significant results year by year.

# Percent change from 2020 to 2024 for both election result difference (pct points) and vaccine rates.
election_pct_change <- vaccine_election_2016 %>% 
  inner_join(vaccine_election_2024, by=c("cntyname" = "cntyname")) %>%
  mutate(pct_change_election_points = 100 * (per_point_diff.y - per_point_diff.x)/per_point_diff.x) %>% 
  select(cntyname, Y_2016, Y_2024, pct_change_election_points)

# Correlation between percent change in infant vaccine rates and percent change in election difference.

# Join both percent change datasets together.
cc_baby_election <- baby_vax_change %>% inner_join(
  election_pct_change, by=c("cntyname" = "cntyname")
) %>% 
  select(cntyname, pct_change_14_24, pct_change_election_points)

# Coefficient is -0.07. Little correlation.
cor(cc_baby_election$pct_change_14_24, cc_baby_election$pct_change_election_points)
## [1] -0.07368496
# Same thing for teen vaccine series. Join both percentage change datasets together.
cc_teen_election <- teen_vax_change %>% inner_join(
  election_pct_change, by=c("cntyname" = "cntyname")
) %>% 
  select(cntyname, pct_change_14_24, pct_change_election_points)

# Coefficient is -0.012. Little correlation.
cor(cc_teen_election$pct_change_14_24, cc_teen_election$pct_change_election_points)
## [1] -0.01242924

MMR vaccine analysis from report card data (only 2024)

# Let's look at MMR vaccination rates as of Dec. 31, 2024 (end of 2024 Q4)
report_card_data %>% 
  filter(
    year == "2024",
    quarter == "Q4"
  ) %>% 
  select(year, quarter, cntyname, covgmmr) %>% 
  arrange(covgmmr) %>% 
  slice(1:10)
##    year quarter  cntyname covgmmr
## 1  2024      Q4    Oscoda    52.8
## 2  2024      Q4  Keweenaw    56.5
## 3  2024      Q4   Gladwin    62.9
## 4  2024      Q4  Mackinac    66.9
## 5  2024      Q4   Sanilac    67.6
## 6  2024      Q4     Clare    67.9
## 7  2024      Q4   Osceola    69.2
## 8  2024      Q4  Houghton    69.3
## 9  2024      Q4 Hillsdale    70.1
## 10 2024      Q4    Branch    70.3
report_card_data %>% 
  filter(
    year == "2024",
    quarter == "Q3"
  ) %>% 
  select(year, quarter, cntyname, covgmmr) %>% 
  arrange(covgmmr) %>% 
  slice(1:10)
##    year quarter cntyname covgmmr
## 1  2024      Q3   Oscoda    56.4
## 2  2024      Q3  Gladwin    61.5
## 3  2024      Q3 Keweenaw    61.9
## 4  2024      Q3  Sanilac    66.8
## 5  2024      Q3 Mackinac    67.5
## 6  2024      Q3    Clare    69.2
## 7  2024      Q3 Houghton    69.2
## 8  2024      Q3     Lake    69.7
## 9  2024      Q3  Osceola    70.6
## 10 2024      Q3   Branch    71.0
report_card_data %>% 
  filter(
    year == "2024",
    quarter == "Q2"
  ) %>% 
  select(year, quarter, cntyname, covgmmr) %>% 
  arrange(covgmmr) %>% 
  slice(1:10)
##    year quarter cntyname covgmmr
## 1  2024      Q2 Keweenaw    57.8
## 2  2024      Q2  Gladwin    64.5
## 3  2024      Q2  Sanilac    66.9
## 4  2024      Q2 Mackinac    67.8
## 5  2024      Q2 Houghton    68.4
## 6  2024      Q2  Osceola    69.0
## 7  2024      Q2     Lake    70.6
## 8  2024      Q2    Clare    70.7
## 9  2024      Q2  Newaygo    71.4
## 10 2024      Q2    Iosco    72.2
# Aggregating MMR vaccines by county 


# Oscoda County had the lowest MMR vaccination rates, followed by Keweenaw, Gladwin, Mackinac and Salinac. None of these counties saw measles cases in 2025.

Potential sources:

Thrishika Bala, epidemiologist for the state of Michigan

Lynn Sutfin, PIO for Michigan DHHS, SutfinL1@michigan.gov

Ryan Malosh at MaloshR@michigan.gov

Data Details

This website explains the different Series numbers for the sets of vaccines: https://www.michigan.gov/mdhhs/adult-child-serv/childrenfamilies/immunizations/data-statistics/countyimmsreportcards

Source: https://www.michigan.gov/mdhhs/adult-child-serv/childrenfamilies/immunizations/data-statistics/countyimmsreportcards

Documentation: https://www.michigan.gov/mdhhs/-/media/Project/Websites/mdhhs/Adult-and-Childrens-Services/Children-and-Families/Immunization-Information/LHD/Immunization-Report-Cards/Report_Card_Defs__FAQs_2021.pdf?rev=93e4ac10b66c4c10abcd3c02a6b55982&hash=EA6122E84CB516E9E6D335CC73BDB5CF